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Abstract 

We study symmetric sleepy random walkers, a model exhibiting an absorbing-state phase 
transition in the conserved directed percolation (CDP) universality class. Unlike most ex- 
amples of this class studied previously, this model possesses a continuously variable control 
parameter, facilitating analysis of critical properties. We study the model using two com- 
plementary approaches: analysis of the numerically exact quasistationary (QS) probability 
distribution on rings of up to 22 sites, and Monte Carlo simulation of systems of up to 32000 
sites. The resulting estimates for critical exponents /3, and z, and the moment ratio 

"1211 = {p 2 )/{p} 2 (p is the activity density), based on finite-size scaling at the critical point, 
are in agreement with previous results for the CDP universality class. We find, however, 
that the approach to the QS regime is characterized by a different value of the dynamic 
exponent z than found in the QS regime. 
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I. INTRODUCTION 



Over the last several decades, phase transitions between an active and an absorb- 
ing state have attracted great interest in statistical physics and related fields jl-4|. 
More recently, experiments on such transitions have been performed 5KZ|. As in 
equilibrium, continuous phase transitions to an absorbing state can be grouped into 
universality classes jsl, 0] . Two classes that have received much attention are directed 
percolation (DP) and conserved directed percolation (CDP), exemplified, respectively, 
by the contact process and the stochastic conserved sandpile (conserved Manna 



model) 10, While the former class is well characterized, and there is a clear, con- 



sistent picture of the scaling behavior, the critical exponents of CDP have not been 
determined to high precision, and there are suggestions of violations of scaling. Thus 
it is of interest to study further examples of this class, and to apply new methods of 
analysis to such models. Absorbing-state transitions have been studied via mean-field 
theory, series expansion {2], renormalization group jsj, perturbation theory 12] and 
numerical simulation. Recently, an analysis based on the exact (numerical) determi- 
nation of the quasistationary (QS) probability distribution was proposed and applied 
to models in the DP class 13|. 

In this paper we study sleepy random walkers (SRW), a Markov process defined on a 
lattice, belonging to the CDP universality class, using exact (numerical) analysis of the 
quasistationary (QS) probability distribution and Monte Carlo simulation. The former 
approach furnishes quite accurate predictions for the contact process; a preliminary 
application to a model in the CDP class yielded less encouraging results, due in part 
to the small system sizes accessible 13]. The smaller number of configurations (for a 
given lattice size and particle density) in the SRW model allows us to study somewhat 
larger systems, leading to improved results in the QS analysis. We study the model 
in extensive Monte Carlo simulations as well, in efforts to better characterize CDP 
critical behavior. A closely related model, activated random walkers (ARW), was 
introduced in |14||; in this case there is no restriction on the number of walkers per 
site. In |14j the principal emphasis was on asymmetric ARW (hopping in one direction 
only); some preliminary evidence for CDP-like behavior of the symmetric version was 
also reported. 



The balance of this paper is organized as follows. In Sec. II we define the model 
and its behavior in mean-field theory. In Sec. Ill we describe how exact QS analysis 
is applied to the model and present the associated results on critical behavior. We 
report our simulation results in Sec. IV, and in Sec. V we present a summary of our 
findings. 

II. MODEL 

The SRW model is defined on a d- dimensional lattice of L d sites with periodic 
boundary conditions. Each site % of the lattice may be in one of three states: empty 
(<7j = 0), occupied by an active particle (a* = 1), or by an inactive particle (<7j = — 1). 
Multiple occupancy is forbidden. Active particles attempt to hop, at unit rate, to 
a nearest-neighbor site. In a hopping move the target site is chosen with uniform 
probability on the set of nearest neighbors, and the move is accepted if and only if the 
target site is vacant. Transitions from o~i = I (active) to a, = — 1 (inactive) occur at a 
rate of A, called the sleeping rate, independent of the states of the other sites. Inactive 
particles cannot hop. A transition from <ji = —1 to Oi = 1 occurs when an active 
neighbor attempts to jump to site i. In this case, the particle that attempted to hop 
returns to its original site, but the particle that was sleeping is activated. Evidently 
this Markovian dynamics conserves the number of particles. Here we consider initial 
configurations in which N particles (all active) are distributed randomly amongst the 
sites, respecting the prohibition of multiple occupancy. (In the ARW model [14] the 
number of particles per site is unrestricted but only an isolated particle can sleep.) 

Let N a denote the number of active particles; any configuration with iV a = is 
absorbing. Thus we define the order parameter as p = N a /N, the fraction of active 
particles. There are two control parameters, the sleeping rate A and the particle 
density ( = N/L d . For £ < 1, the particle number is a nontrivial conserved quantity, 
and we expect the model to belong to the CDP universality class. (For N = L particle 
conservation follows trivially from the conservation of site number, and the model is 
equivalent to the contact process, belonging to the DP class.) 

An advantage of this model is the presence of a continuously variable control pa- 
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rameter, A. In the stochastic sandpile 10|, LL3J, the control parameter, (, cannot only 
be varied in increments of 1/L d , which tends to complicate the determination of crit- 
ical properties. (A given value of ( is accessible only for a restricted set of system 
sizes.) We shall therefore fix the particle density and vary A to locate the critical point. 
Since the QS distribution analysis depends on applying finite-size scaling analysis, we 
use a £ = 1/2, which is accessible in all systems with L even. 

Mean field (MF) analysis yields the following equation of motion for the fraction 
of active particles: 

which is analogous to the MF equation for the contact process (CP) |2| if we identify 
£ and A with the creation and annihilation rates, respectively, in the CP. One sees 
immediately that at this level of approximation, an active stationary state exists only 
for A < A c = C, in which case the stationary order parameter is p = £ — A. Although 
the MF analysis is certainly not reliable in detail, it is reasonable to expect that the 
model exhibits a continuous phase transition, and that A c is an increasing function of 

c 



III. QUASISTATIONARY ANALYSIS 

A. Quasistationary probability distribution 

In [3], one of us proposed a method for studying absorbing-state phase transitions 
based on numerical determination of the quasistationary probability distribution, that 
is, the asymptotic distribution, conditioned on survival. With the essentially exact 
QS properties in hand, one may apply finite-size scaling analysis to estimate critical 
properties. 

Let p c = \im t ^oo Pc(t)/P(t) denote the QS probability of configuration c, where 
p c (t) is the probability at time t and P(t) is the survival probability, i.e., that proba- 
bility that the absorbing state has not been visited up to time t. The QS distribution 
is normalized so: Yl c Pc = 1> where the sum is over nonabsorbing configurations only 
(the QS probability of any absorbing configuration is zero by definition). Given the 
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set of all configurations (including absorbing ones) and the values of all transition 
rates to c /,c (from c to c'), we construct the QS distribution via the iterative scheme 
demonstrated in [151 ]: 



p' c = ap c + (1 - a) — (2) 

w c -r a 

Here r c = ^2 C , w C)C 'p c ' is the probability flux (in the master equation) into state 
c, r a is the flux to the absorbing state (l/r a gives the lifetime of QS state), and 
w c = Ylc' w c',c is the total rate of transitions out of state c. The parameter a can 
take any value between and 1 (in practice we use a = 0.1). Following each iteration, 
the resulting distribution p' c is normalized by multiplying each probability by / = 
VEcPcl 1 Starting from an initial guess (for example, a distribution uniform on 
the set of nonabsorbing configurations), this scheme rapidly converges to the QS 
distribution. 

Since the number of configurations and transitions grows very rapidly with system 
size, we use a computational algorithm for their enumeration. To begin, we enumer- 
ate all configurations of L/2 particles on a ring of L sites (recall that each particle 
must occupy a distinct site). Configurations that differ only by a lattice translation 
or reflection are treated as equivalent. Thus the space of configurations is divided 
into equivalence classes C. For each class we store one representative configuration, 
and the number |C| of configurations in the class, which we call its weight. Each 
configuration is determined by (1) the particle positions and (2) the state (active or 
sleeping) of each particle. If we ignore the particle states, the particle positions define 
the basic configuration; each basic configuration corresponds to a series of configura- 
tions c. One such configuration has all particles active, while others have 1, 2, L/2 
inactive particles; the one with all particles inactive is absorbing. Once the set of 
basic classes has been enumerated, we enumerate the classes with n p = 0, 1, ...,L/2, 
inactive particles, and their associated weights. 

Next, we enumerate all transitions between configurations. We visit each equiva- 
lence class C in turn, and enumerate all the manners in which C arises in a transition 
(due to particle hopping, inactivation, or activation) from an antecedent configuration 
in some class C. Each transition is characterized by a rate u>c,c an d by an associated 
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weight, mc t c'- (The latter is needed because in certain cases, two or more distinct 
transitions to the same class C have antecedent configurations belonging to the same 
class, C.) Given the set of classes and transitions, and associated rates and weights, 
we can iterate the relation given above to determine the QS probability distribution on 
the set of nonabsorbing classes. (The sums are now over classes, with normalization 
taking the form Y2c \C\Pc = !•) 

We determine the QS distribution on rings of size L — 6, 8, 10, 22. For L — 22 
the total number of equivalence classes is N con f = 32 842 718, and the number of 
transitions involving hopping and sleeping are = 265 512 131 and N s = 180 594 624, 
respectively. Our criterion for convergence of Eq. (J2J) is that the sum of all absolute 
differences between the probabilities pc and p' c at successive iterations be smaller than 

io- 15 . 



B. Critical properties 



Extracting estimates for critical properties from results for small systems depends 
on finite-size scaling (FSS) analysis 16|, Il7j . The FSS hypothesis implies that the 
order parameter follows p(A,L) oc L~^^ U± TZ(L 1 ^ U± A), where A = (A c — A)/A c and 1Z 
is a scaling function. (Note that in the SRW model the active phase corresponds to 
A < A c .) The QS order parameter is given by p = (L/2) -1 J2 c p c N a!C , with iV 0>c the 
number of active particles in configuration c. To find the critical exponent (3 = 0/v±, 
we seek crossings of the quantities [18] , 



S L (X) = 



ln[p(A,L + l)/p(A,L-l)] 
ln[(L+l)/(L-l)] 



(3) 



for successive pairs of system sizes. Let Sl + i(A) = Sl_i(A) = f3(L) for A = \s,l- The 
crossing values \s,l and (3(L) are expected to converge to A c and (3, respectively, as 
L — t- oo. 

To estimate the dynamic exponent z, we determine the QS probability flux to the 
absorbing state (i.e., the inverse lifetime), which follows r a oc L^ z J r (AL 1 ^ u ), with T 
another scaling function. The crossings of 
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p m _ ln[r a (A,L-l)/r a (A,L + l)] 
L( } = M(L + 1)/(L-1)] (4) 



furnish a series of estimates, Zl- As in the case of Sl above, the A values, \r,l, at the 
crossings are expected to converge to Ac- 
Critical behavior at an absor bing state phase transition is also characterized by 
order-parameter moment ratios [19]]. Let denote the k-th moment of the or- 
der parameter. The scaling property of the QS probability distribution leads to the 
asymptotic size- invar iance of moment ratios of the form m n j '(m*m^) for ir + js = n, 
at the critical point. (Although not, strictly speaking, a ratio, the product m_im of 
the first positive and negative moments follows the same general scheme.) We analyze 
the ratios 771211 = 7712/m 2 , TO3111 = m^/mf, m_im, and the reduced fourth cumulant, 
or kurtosis q. The latter is defined so: q = K^/K^, where K 2 = m 2 — m\ = var(p) 
and K± = — Am^uii — 3m 2 , + V2m<ivn\ — 6mf. The A values marking the crossings of 
the moment ratios (for system sizes L and L + 2) are once again expected to converge 
to A c as L — > 00. The values of the moment ratios and q at the critical point are uni- 
versal quantities, determined by the scaling form of the order-parameter probability 
distribution 
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20], and so are useful in identifying the universality class. 
A further quantity of interest is the scaled variance of the order parameter x — 
L d (p 2 — p 2 ) which is expected to diverge as |A — A c |~ 7 . (In equilibrium systems, x is 
proportional to the susceptibility). In a system of size L, x exhibits a maximum at a 
sleeping rate we denote X x ,l- FSS predicts that at the critical point x ^ L J / U± . 

While the quantities mentioned above furnish the exponent ratios (3/v±, v\\/v±. = -2, 
and 7 = r y/v±, it is also possible to estimate i>± directly. FSS implies that m-in — 
1Z(/\L l l u ^), where 1Z is a scaling function. Thus r' = \dm2n,L/d\\\ c oc L l l u ^. The 
derivatives of other moment ratios, of In p, and of In r a scale in an analogous manner. 

Given a series of estimates for the critical point (or for a critical exponent, or a 
moment ratio), associated with a sequence of sizes L, we extrapolate to infinite size 



using polynomial fits and the Bulirsch-Stoer (BST) procedure 24|. In the contact 



process [13], quantities such as Xs,l and f3(L) vary quite systematically with system 



size, leading to precise estimates for critical values via BST extrapolation. 
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Our first task is to determine the critical sleeping rate A c ; to this end we analyze 
the crossings of S, R, and the moment ratios. In a preliminary analysis the crossings 
are determined graphically; the general tendencies are shown in Fig. [TJ Subsequently, 
we refine these estimates by calculating these quantities at intervals of 5 A = 10 -5 , 
at 20 points around each estimated crossing, and determine the crossing values to 



a precision of 10 -12 or better using Neville's algorithm 23|]. Figure [2] shows the 
results for crossings of S^. This quantity is unusual in that the crossing values are 
nonmonotonic; for the other quantities studied, the L-dependence is monotonic over 
the accessible range of system sizes (see Fig. Hj). 

In BST extrapolation [24j the limiting value of a sequence Tj (J = 1,2,3, ...), 
is estimated on the basis of the first N terms via the recurrence relations 



j>(n) _ rp(n+l) _|_ /rp(n+l) _ rp 
m m—1 ' \ m—1 i 



m—1/ 



h, 



L n+m 



rp(n+l) rp(n) 

m—1 m—1 



rp(n+l) rp 

1 m-l 1 - 



(n+1) 
2 



-i -1 



(5) 



where, for j = 1, N, T_l = 0, Tr) 3 ' = Tj, and hj is a sequence converging to zero 
as j oo. (Here hj = 1/Lj, where Lj is the system size, or, for crossings, the mean 
value of the two system sizes involved.) 

The BST procedure includes a free parameter, u, which can be adjusted to improve 
convergence. We use a convergence criterion similar to one employed in analyses of 
series expansions via Pade approximants 25|, |26[ , in which, varying some parameter, 
one seeks concordance amongst the estimates furnished by various approximants. In 
the present case, given iV values Tj, each associated with inverse system size hj, we 
calculate N + 1 estimates, one using the full set, and iV others obtained by removing 
one point, (hk,T^), from the set. We search for values of u that minimize the differ- 
ences between the various estimates. Sweeping the interval [0,5], we find that each 
quantity studied (/3, z, moment ratios, and the associated estimates for A c ) exhibits 
one or more crossings at which all iV + 1 estimates are equal to within numerical 
precision. Figure [3], for \ c $ (the estimate for A c derived from the crossing values 
\s,l)i illustrates the typical behavior. In this case there are four crossings, which fall 
at u = 1.051380, 1.658703, 2.109962, and 2.550852; the associated values of A c , 5 are 



S 




FIG. 1: QS analysis: R, S, m.211, w-3111, m_im in the neighborhood the crossings. The 
insets show these quantities over a larger range of A values. The lower-right panel shows 
the crossing values for m2n,m_imi and 7713111 (lower to upper) along with the extrapolated 
(L — > 00) values. 

0.0904577, 0.0903878, 0.0902354, and 0.0900773. 

To choose among the values when there are multiple crossings, we note that u in 
the BST procedure is effectively a correction to scaling exponent. An independent 
estimate for this exponent can be obtained via a least-squares fit to the data using a 
double power-law form, for example, 

A B 

X s ,l = X c ,s + J^ + J^ (6) 
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0.245 




0.00 0.02 0.04 0.06 0.08 0.10 0.12 
1/L 

FIG. 2: Crossing values (3{L) versus 1/L. The inset is a similar plot of the crossing points 
Xs,L- The leftmost points are the extrapolated values using BST; curves are splines to the 
data and the BST extrapolations, intended as a guide to the eye. Error bars are smaller 
than the symbols. 

with jj2 > Hi- (The best-fit parameters A, B, y\ and yi are determined by minimizing 
the variance of the differences 5l = ^s,l ~ AL~ Vl — BL~ m .) This yields y\ = 2.02, 
leading us to take the average of the two values associated with the BST crossings 
nearest yi, resulting in A C) g = 0.09016. A similar procedure is used to obtain the other 
estimates listed in Tables fl] and flU (We note that the apparent correction to scaling 
exponent y\ falls in the range 2.02-2.23 for the crossing values of A, and in the range 
1.1-1.6 for the associated quantities (3, z, and the moment ratios.) 

Tables [J and [TT] include polynomial extrapolations as alternative estimates for the 
quantities of interest. (In this case we fit the data to a polynomial in 1/L, using the 
highest possible degree. Polynomials of degree one or two smaller than maximum yield 
very similar results.) We adopt the mean of the BST and polynomial extrapolations 
as our best estimate, and adopt the difference between the two results as a rough 
estimate of the associated uncertainty. The situation is particularly favorable for 
determining A c since we have five independent estimates. The average of the BST 
results is A c = 0.08996(7) while that from the polynomial fits is 0.09007(10), leading 
to our best estimate of A c = 0.09002(10). (Figures in parentheses denote uncertainties, 



10 



given as one standard deviation.) The estimates for A c and the moment ratios are 
consistent with simulation results, whereas those for (3/i/± and z are not. (A detailed 
comparison is given in the following section.) Figure H] shows the finite-size data and 
BST extrapolations for the various crossing values, A M £. 




FIG. 3: BST estimates (see text) for A C) 5 versus extrapolation parameter u. Main graph: 
detail of interval [2,3]; inset: the full interval of study. 
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FIG. 4: Crossing values A. 5 l for Rl, Sl, w-211, ""13111 and m-\m\ (upper to lower). 

The values for the kurtosis at A c approach a limit of q c = —0.454(2) (see Fig. 
For a given system size, q(X, L) exhibits a minimum in the vicinity of A c , as also 
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observed in the contact process 131 ] . Extrapolating q(X,L) to L — > oo for a series of 
values near A c , we obtain a function that exhibits a minimum near A = 0.0947 (the 
minimum value is -0.521). This implies that the minimum falls near, but not at the 
critical point, a conclusion supported by the simulation data reported in the following 
section. Turning to the scaled order-parameter variance \i we observe pronounced 
maxima even in small systems, as illustrated in Fig. El Estimates for 7 obtained from 
a local-slopes analysis of x a t the critical point A c are listed in Table [TT1 (The local 
slope is defined so: j{L) = \n[ X (L + l)/x(L - l)]/ln[(L + \)/{L - 1)].) 




0.060 0.075 0.090 0.105 0.120 0.135 0.150 



X 



FIG. 5: Kurtosis q vs A for sizes 6, 8, 22 (lower to upper). Inset: values for q{L) at the 
critical point and our estimate for q\ Ct00 . 



To estimate u±, we determine r' = \dm2u,L/ by constructing linear fits to the 
data on the interval 0.08992 < A < 0.09012, using an increment of AA = 10~ 5 . Since 
the graph of In r' versus In L shows significant curvature, we analyze the local slopes, 
z/j_(L) = ln[(L — 1)/(L + 1)]/ ln(r' L _ 1 /r' L+1 ). BST extrapolation of the latter yields 
u± = 1.293(5). We obtain independent estimates for u± using the order parameter and 
flux of probability to the absorbing state, r a , as described above, yielding 1.285(10) 
and 1.2644(2), respectively. On the basis of these results, we estimate u± = 1.28(1). 
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FIG. 6: Scaled order parameter variance \ versus A for sizes 6, 8, 22 (lower to upper). 
Inset: Estimates for 7/^j_ obtained via local-slopes analysis, and the extrapolated (infinite- 
size) value. 



TABLE I: Estimates for the critical sleeping rate A c obtained via analysis of the QS proba- 
bility distribution using BST extrapolation and polynomial fits. 



Quantity 


BST 


Polynomial 


A C ,S 


0.09016 


0.09039 


Ac,/? 


0.08973 


0.08993 


A c ,211 


0.08999 


0.09008 


A c ,3111 


0.08995 


0.09016 


A c ,-n 


0.08998 


0.08979 



IV. MONTE CARLO SIMULATIONS 



A. Simulation methods 



We perform extensive simulations of the SRW model using both conventional and 
quasistationary (QS) methods. Quasistationary simulations |2l|, |22| have proven to 
be an efficient method for studying absorbing state phase transitions, allowing one to 
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TABLE II: Estimates for critical properties obtained via analysis of the QS probability 
distribution using BST extrapolation and polynomial fits. 



(Quantity 


tib 1 


T"> 1 "1 

Polynomial 


Best Est. 


P/v± 


0.2418 


0.2405 


0.241(1) 


z 


1.668 


1.660 


1.664(4) 


y/v ± 


0.5428(1) 


0.53942(1) 


0.541(2) 


"1211 


1.1412 


1.1422 


1.1417(5) 


"2-3111 


1.4151 


1.4249 


1.420(5) 




1.2965 


1.3203 


1.308(12) 


Q 


-0.460(5) 


-0.454(4) 


-0.457(6) 



obtain results of a given precision with an order of magnitude less CPU time than 
in conventional simulations. The method samples the QS distribution defined in the 
preceding section using a list of configurations saved during the evolution; when a visit 
to the absorbing state is imminent, the system is instead placed in a configuration 
chosen at random from the list. A detailed explanation of the method is given in 21|. 

We perform QS simulations using system sizes L = 100, 200, 400,..., 32000. Each 
realization of the process runs for T = 10 9 time units, with the first Tr = 10 8 time 
units discarded to ensure all transients have been eliminated. (Our time unit is 
defined below.) We use 1000 saved configurations; the replacement probability (i.e., 
for replacing one of the configurations on the list with the current one) is p rep = 10 -5 . 
Our choice of p rep is guided by the condition T > t m > r, where t m = M/p rep is 
the mean time that a configuration remains on the list and r is the mean lifetime 
in the QS state. The latter is estimated as r = (T — T R )/N abs , where N abs is the 
number of (attempted) visits to the absorbing state for t > Tr. During the initial 
relaxation period (t < Tr) we use p rep = 10~ 2 to eliminate the memory of the initial 
configuration. For each value of A studied, we calculate the mean and statistical 
uncertainties (given as one standard deviation) over Nr = 20 independent realizations; 
for L = 16000 we use N R = 40. 

At each step of the simulation we select the particle involved from a list of active 
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particles. The time increment associated with each step is At = 1/N a , where N a is 
the number of active particles just prior to the event. For t > Tr, we accumulate 
a histogram h(N a ) of the time during which there are exactly N a active particles. 
The normalized histogram is our best estimate for the probability distribution P(N a ), 
from which we may determine any desired moment of the order parameter p = N a /N . 
The QS lifetime r may also be obtained from P{N a ) via the relation = 1/[AP(1)], 
where the subscript h serves only to distinguish this from the value r found using the 
mean time between visits to the absorbing state. 

B. Scaling at the critical point 

To determine A c , we first locate the crossings of the moment ratios m-iw (defined 
in Sec. Ill), for pairs of consecutive system sizes, and obtain a preliminary estimate 
by extrapolating the crossing values to L — > oo. We then study larger systems using 
A values close to our preliminary estimate. Using these results, we determine the 
critical value via the familiar finite-size scaling criteria p ~ L~P /vx and r ~ L z ', and 
the condition that 777.211 approach a finite limiting value as L increases. In logarithmic 
plots, p and r exhibit upward (downward) curvature for A < A c (> A c ) as illustrated 
in Figs. [7J and [HJ off-critical values are also readily identified in plots of m 2 n versus 
1/L. Using these criteria we obtain A c = 0.090085(12). Of note are the strong finite- 
size corrections (evident in the insets of Figs. [7] and El), for L < 1000. Indeed, our 
estimates for critical exponents and moment ratios are obtained using only the data 
for L > 1000. 

We turn now to FSS estimates of critical exponents. Using the data for L > 1000 
we obtain /3/Vj_ = 0.212(6) from analysis of p, and 0.217(10) from analysis of m_i, 
which as noted in Sec. Ill, is expected to follow m^i(L, A c ) oc L^' u± . Analysis of the 
QS lifetime using r and yields z = 1.50(4) and z = 1.51(4) respectively, while the 
data for the moment ratio yield m^w^ = 1.141(8). Restricting the analysis to the 
data for L > 2000, or for L > 4000, yields estimates for /3/v±, 771211,0 an d z consistent 
with the values cited above, but with somewhat larger uncertainties. Analysis of x 
yields 7/z/_i_ = 0.58(1). In all cases the chief contribution to the uncertainty is due to 
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FIG. 7: Order parameter p versus system size L. Inset: L°' 2l2 p versus L. Lines are a guide 
to the eye. 
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FIG. 8: Lifetime versus system size L. Inset: L versus L. Lines are a guide to the 

eye. 

the uncertainty in A c itself. 

To estimate the exponent u± directly, we apply the method used in Sec. IHIt 
calculating the derivatives r' of quantities such as In p and In r with respect to A near 
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the critical point. Using simulation data for A = 0.090073, 0.090085, and 0.090097, we 
construct a linear fit to estimate r' at A c . In Fig [9] we plot the values for r' obtained 
via analysis of \d\np/d\\, dm,2ii/d\ and |dlnr/dA|. The estimates obtained using 
the data for L > 1000 are listed in Table IHIt based on these results we estimate 
u± = 1.31(4). Since the values obtained using different quantities are quite different, 
this exponent is not determined to good precision. 



T3 

o 

"D 




1000 10000 



L 

FIG. 9: Derivatives r' = \dQ/d\\\ c for m 2 n (squares), In/? (circles) and lnr (triangle). 
Lines are linear fits to the data, with slopes of 0.746(37), 0.749(16), and 0.80(3) (lower to 
upper). 



TABLE III: Estimates for i>± obtained from the derivatives of m2n, hip and lnr. 



r> = \dQ/d\\ Xc 


Q = m 2 H 


Q = Hp) 


Q = ln(r) 


1000 < L < 32000 


1.34(7) 


1.34(4) 


1.25(5) 



Next we examine the QS probability distribution Pl{p, A) for the fraction of active 
particles, p = N a /N. At the critical point, this distribution is expected to take the 
scaling form j^], 

Pl(p,X c ) = j^ I P(p/(p)), (7) 
where P is a normalized scaling function. (The prefactor arises from normalization of 
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Pl, with N = L/2.) Figure [TD] is a scaling plot of the QS probability distribution at 
the critical point, showing evidence of a data collapse, albeit with significant finite- 
size corrections (the maximal scatter of the values is about 2%); the collapse is quite 
good for the two largest system sizes. Also shown are scaled probability distributions 
for the one-dimensional conserved restricted stochastic sandpile j^], showing good 
overall agreement. 

Figure [TT] shows the behavior of the order-parameter moment ratios m 2 n, m 3 m, 
the reduced fourth cumulant q, and m_im, defined in Sec. III. Our estimates for the 
critical values are given in Table PVT) . 




p/<p> 

FIG. 10: Main figure: scaling function P = ((p)L/2)Px(p/(p)) versus p/(p) for system sizes 
L = 1000, 2000,. ..,L = 32 000 (lower to upper). The dotted and dashed curves show the 
corresponding result for the one-dimensional restricted stochastic sandpile, for L = 20 000 
and 50000, respectively. Lower inset: detail of the region in which P takes its maximum. 
Upper inset: unsealed data (system sizes increasing from upper to lower). 

C. Off-critical scaling behavior 

It is of interest to study the scaling of the order parameter, of the scaled variance 
and of the lifetime, away from the critical point. Although off-critical scaling proper- 
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FIG. 11: Moment ratios (a) 771211; (b) mzxn', (c) reduced fourth cumulant q; and (d) m_im 
versus system size L. 



ties have been amply verified 
such as the contact process 



'or models in the directed percolation universality class, 



271 ]. finite-size scaling and associated data collapse 



of the order parameter is more problematic in one- dimensional stochastic sandpile 



models 
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FSS analysis implies that the order parameter take the form p(A, L) ~ 
Ij-P/"± J(AL 1 ' i/ - l ), where the scaling function is f(x) oc x 13 for x 1. In the in- 
active phase, A > A c , the number of active particles in the QS regime is 0(1), so that 



p oc L 1 , leading to /(x) 



0-u ± 



. Alternatively, writing p(A, L) ~ A /3 /i(AL 1 /^ 



the scaling function must satisfy /i(x) oc x 13 for a; — > 0, while in the inactive 
phase h(x) oc |x| _i/± . FSS analysis predicts that the scaled variance take the form 
X(L,A) ~ L rf/u - L g(AL 1/u - L ), with g(x) oc x" 7 for x > 1 and #(x) oc |x|" (i/±+7) in the 
inactive phase. 

Figure [12] shows a good data collapse of the order parameter in the forms p* = 



]_,Pl v ±p and p = pA - ^, and of the scaled order-parameter variance x* = xL~~ / ^ u± , as 
functions of A* = AL 1 ^ U± , using data for system sizes from L = 100 to 32000. The 
exponents associated with the best collapses are listed in Table [TV] Based on these 
results, we estimate f}/u± = 0.216(4), = 0.29(1) and 7/1/j. = 0.57(1). We note that 
the values found for u± in the active and inactive regimes differ slightly, leading to 
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best estimate u± = 1.30(4). 
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FIG. 12: Scaled order parameter p* = pL^/ u± - (upper left) and p = pA~@ (upper right), 
scaled variance x* — X-^~ 7 ^ x (lower left) and the scaled lifetime r* = L~ z t (lower right) 
versus scaled distance from critical point, A* = AL 1//l/J -, in the active and inactive phases 
(upper and lower set of points, respectively). The best-fit exponents associated with the 
data collapses are given in Table HVl The slopes of the dashed and dotted lines represent the 
power laws exhibited by the scaling functions in the active and inactive phase, respectively. 

In the active phase, A < A c , a data collapse of p* versus A* is obtained over about 
four orders of magnitude in A*. Interestingly, such a data collapse is only observed 
over a much smaller interval - about one order of magnitude - in stochastic sandpiles for 
a comparable range of lattice sizes 29_|, |30j. A linear fit to the data (using all sizes) for 



A* > 10 yields (3 = 0.293(2). Using p versus A* we found (3 = 0.288(3), including the 
data for A* < 0.2. The data for x* (for A* > 10) yield 7 = 0.733(7). The power laws 
associated with these exponents are represented in Figure [12] by dashed lines. Using 
the hyperscaling relation 7 = duj_ — 2(3, the latter results and the exponents used in 
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the collapses in the active phase, one finds (a) (3 = (1.32(3) — 0.733(7))/2 = 0.29(2); 
(b) p/ v± = (1 - t/i/±)/2 = (1 - 0.568)/2 = 0.216(5); and (c) u± = 7 + 2(3 = 
0.733(7) + 0.580(8) = 1.31(2). These predictions are consistent with the values used 
in the collapses and with those obtained at the critical point. 

TABLE IV: Critical exponent estimates from off-critical simulations. For each scaling rela- 
tion (first column), we list the associated exponent obtained via a fit to the data (second 
column). The third and fourth columns give the exponents obtained via data collapse, x 
denotes the argument of the relevant scaling function. 



Active phase 


p* oc 


P = 0.293(2) 


P/u± = 0.218(5) 


v ± = 1.34(3) 


p oc x~P 


P = 0.288(3) 


P = 0.293(3) 


v ± = 1.32(3) 


X* oc x~ 7 


7 = 0.733(7) 


= 0.568(8) 


v± = 1.32(3) 


Inactive phase 


p * „ x P-v± 


P - v L = -0.95(4) 


P/u± = 0.214(2) 


v ± = 1.26(2) 


p ~ x~@ 


P = 0.29(1) 


u ± = 1.22(5) 


v ± = 1.26(3) 


X* ~ x~^ ±+7 ) 


v± +7 = 1.85(4) 


l/v± = 0.57(1) 


v ± = 1.28(3) 


* — Vw 
T ~ X " 


= 1.86(6) 


v\\/v ± = 1.53(3) 


v ± = 1.30(2) 



Similarly, the scaled quantities p*, p and x* exhibit a good collapse in the inactive 
phase, for all system sizes studied, as shown in Fig. [12j Using p*, a linear fit to the 
data for |A*| > 20 yields v± — (3 = —0.95(4). Using the result from the best collapse 
we have = 0.95(4) — 1.26(2) = 0.29(4). Analyzing the data collapse for p we find 
P = 0.29(1) in the small- 1 A* | regime. In the opposite limit we obtain u± = 1.22(5) 
as illustrated in the dotted lines in Fig. [T5J Finally, the data collapse of x* versus 
A* leads to u± + 7 = 1.85(4). Combining these results, the hyperscaling relation 
used above, and the values from the best collapses shown in Table HVl for x* 1 one 
can predict: (a) j/u ± = 0.57(1) 7 = 0.57(1) x 1.28(3) = 0.73(3), (b) p/u± = 
(1 - 7/i/±)/2 = (1 - 0.57(l))/2 = 0.215(5) and (c) 7 = 1.85(4) - 1.28(3) = 0.57(5). 
The first two relations are consistent with the exponent values obtained previously, 
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while the third conflicts with the value for 7 found for A < A c . This may reflect a 
violation of scaling, but the possibility that our study does not probe sufficiently deep 
into the inactive regime, for which one expects x ~ cannot be discarded. 

In the inactive phase, the lifetime r is expected to follow r(A, L) = 
|A| _I/ ii(jr(AL 1 / I '- L ), with the scaling function G(x) oc implying a data collapse if 

we plot r* = L~ z t versus A*. As shown in Fig. [T2J the collapse here is much poorer 
than in the other cases. A linear fit to the data for L < 8000 and A* > 10 furnishes 
v\\ = 1.86(6). This result is in concordance with the scaling relation z = v\\/v±, if 
we employ the values found at the critical point (z = 1.50(4) and u± = 1.34(6)) and 
the best-collapse values (z = 1.53(3) and u±_ = 1.30(2)); the latter yield the value 
v\\ = 1.53(3) x 1.30(2) = 1.99(8), consistent (to within uncertainty) with the value 
found via collapse. 

D. Approach to the quasistationary regime 

Yet another aspect of scaling at an absorbing-state phase transition involves the 
approach to the QS regime, starting from a maximally active initial condition. The 
quantities of interest are the time- dependent activity density p(t) and the moment 
ratio 777,211 (t). At short times (i.e., before the QS regime is attained), at the critical 
point, the expected scaling behavior for the order parameter is p ~ t~ s f(t/r) where 
the scaling function f(x) ~ x s , for i>1, with 5 = /3/v\\. One expects 777.211 — 1 ~ t x l z 
for t < r. 

To probe this regime we perform conventional simulations for the same system 
sizes as used in QS simulations; averages are calculated over Nr = 1000 independent 
realizations; each runs until the system reaches the absorbing state or attains a max- 
imum time, t max = 10 8 . We use the critical point value A c = 0.090085 found in QS 
simulations. We determine p(t) and m 211 (t) as averages over surviving realizations 
on time intervals that for large t represent uniform increments of ln(t), a procedure 
known as logarithmic binning. 

The evolution of p as a function of t is shown in Fig. [13j Unlike the contact process, 
in which p follows a simple power law before attaining the QS regime, in the SRW 
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model the relaxation is more complicated. Absence of simple power-law relaxation of 
the activity density has also been noted for stochastic sandpiles 35j. We calculate 



the local slope S(t) via piece- wise linear fits to the data for hip versus hit on the 
interval [t/a,at], with a = 1.76 (that is, about 30 data points, with an increment 
of 0.1 in lnt). The inset of Fig. [TBI shows 5(t) decreasing systematically with time, 
from around 0.145 to approximately 0.12 over the interval studied. (Note that the 
data for longer times come from larger systems.) If we associate with the short- and 
long-time regimes the values 5 s h or t ~ 0.143(3) and 5i ong ~ 0.121(3), then the scaling 
relation 5 = = ft/(y±_z) yields 5 = 0.212/1.50 = 0.14, in good agreement with 
the short-time value. Plotting p* = L@/ U± p as a function of t* = t/L z we observe a 
good collapse of the data for different system sizes in the first regime using z = 1.50, 
and in the second regime if we use z = 1.72. In both cases the best collapse is obtained 
using p/u ± = 0.215. 




10° 10 1 10 2 10 3 10 4 10 5 10 6 10 7 10 8 

t 



FIG. 13: Active particle density p versus t. System sizes L = 100, 200, 400,... ,L = 32000 
(from upper to lower). The inset shows the local slope of p(t) versus 1/t. 



In contrast to the complicated behavior of the activity density, the quantity 
m 2ii{t) — 1 exhibits simple power-law scaling over four or more decades. Using 
t* = t/L z with z = 1.71, we obtain a good data collapse for all system sizes studied, 
as shown in Fig. HU The relation m%u — 1 ~ t 1 ^ yields z = 1.70(1), consistent with 
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the collapse values for 771211 and p(t) in the long-time regime, but clearly inconsis- 
tent with z = 1.50(4) found using FSS analysis in the QS regime. Curiously, m 2 n 
shows no hint of the crossover exhibited by the order parameter. Note that the value 
5 = 0.121(2) associated with the second regime of the order parameter is consistent 
with /3/v± = 0.212 and z = 1.71. Thus the scaling of m 2 n follows, from the beginning, 
that observed in p in the long-time regime. 
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FIG. 14: Scaling plot of 771211 — 1 versus t* = t/L z at the critical point, using z = 1.70, for 
system sizes up to 32 000. The inset shows the unsealed data for 100 < L < 32000 (from 
top to bottom). The slope of the straight lines is 0.587. 

E. Comparison of exact QS and simulation results 

In Table W\ we compare results obtained via exact analysis of small systems (QSA) 
and simulation. The results are consistent to within uncertainty except for the dynamic 
exponent z. The QSA predictions for j3/v± and (especially) z seem less reliable than 
those derived from simulation. On the other hand, the QSA estimates for moment 
rations are nominally of higher precision than the simulation results. Care must be 
exercised, however, since the QSA analysis may be subject to relatively large finite- 
size corrections. (The present study and previous works on models in the CDP class 
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suggest that corrections to scaling and finite-size effects are stronger for this class than 
for the contact process.) Table IVl also includes results on sandpiles and the conserved 
lattice gas. The agreement between these studies and the present work is quite good, 
leaving little doubt that the SRW model belongs to the same universality class as 
stochastic conserved sandpiles and conserved directed percolation. 



TABLE V: Comparison of critical properties of the SRW model found via exact analysis of small 
systems (QSA) and Monte Carlo simulation (MC), and results from previous studies (Prev) on models 



in the CDP class: a: restricted stochastic sandpile [29(; b: conserved lattice gas 
stochastic sandpile 34 1. 



33[; c: restricted 



Quantity 


QSA 


MC 


Prev 


A c 


0.09002(10) 


0.090085(12) 




P/"± 


0.241(1) 


0.212(6) 


0.213(6) a 


a 




0.290(4) 


0.289(12) a 




0.541(2) 


0.58(1) 


0.55(1) c 


v±_ 


1.28(1) 


1.33(5) 


1.36(2) a 


z 


1.664(4) 


1.50(4) 


1.55(3) b 


m 2 u 


1.142(1) 


1.141(8) 


1.142(8) a 


m 3 m 


1.420(5) 


1.415(26) 


1.425(25) c 


m^im 


1.308(12) 


1.327(27) 


1.332(10) c 


Qc 


-0.454(2) 


-0.47(3) 


-0.46(3) c 



V. DISCUSSION 

We study sleepy random walkers (SRW) in one dimension using (numerically) exact 
quasistationary analysis of small systems and Monte Carlo simulation. Based on 
considerations of symmetry and conserved quantities, one expects the SRW process 
to belong to the conserved directed percolation (CDP) class, typified by conserved 
stochastic sandpiles. Our results for critical exponents and moment ratios support 
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this conclusion. Different from most examples of the CDP class studied until now, 
the SRW process includes a continuously-variable control parameter which facilitates 
simulation and numerical analysis. 

The present work represents a further test of the exact QS analysis proposed in [K 
In addition to locating the critical point with good precision, QS analysis furnishes 
fair results for the critical exponent u± and quite good estimates for the moment ratios 
771211 an d 777.3111- While somewhat better than the preliminary study of a model in 
the CDP class, QSA predictions for the SRW model are not of the quality obtained 
for the contact process 13J. This appears to be connected with the stronger finite- 
size effects and corrections to scaling observed for models in the CDP class. Exact 
analysis of QS properties is nevertheless a useful complement to simulation, as in this 
method the long-time behavior (conditioned on survival) is surely attained, whereas 
simulations are subject to the nagging possibility of insufficient relaxation time. The 
QS calculations typically require (for the largest system studied) several days on a 
reasonably fast computer, that is, a small fraction of the time invested in a simulation 
study. 

In Monte Carlo simulations, we test various scaling relations via data collapse, 
in both the sub- and supercritical regimes, and compare the resulting critical expo- 
nents with those obtained via finite-size scaling at the critical point. The results are 
generally consistent between the regimes, as well as with those of previous studies 
of stochastic conserved sandpiles. Despite the general agreement, we identify several 
inconsistencies and examples of anomalous behavior. First, the estimates for the ex- 
ponent i>± in the inactive phase are significantly smaller (by roughly 5%) than those 
obtained at the critical point or in the active phase. Second, and more significantly, 
the relaxation of the order parameter to its quasistationary value at the critical point 
is marked by two apparent scaling regimes, with associated exponents z s h or t = 1-50, 
8 s hort = 0.143, z long = 1.71, and 5i ong = 0.121, respectively. Paradoxically, z short 
is close to the dynamic exponent characterizing finite-size scaling in the asymptotic 
long-time (i.e., quasistationary) regime. A similar two-regime relaxation is observed 



in the one-dimensional conserved stochastic sandpile 35j; the latter study reports 



z iong = 1-75(3) and somewhat larger values for the exponents 5. Despite these mi- 
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nor numerical differences it seems likely that anomalous relaxation is a characteristic 
of the CDP class in general. Finally, we have noted a possible violation of scaling 
associated with the scaled variance of the order parameter x i n the inactive phase. 

Given the complex pattern of relaxation of the order parameter, it is surprising that 
man — 1, another quantity expected to exhibit power-law scaling at the critical point, 
in fact shows a near-perfect data collapse in a single scaling regime that corresponds 
essentially to the two scaling regimes of the order parameter. The associated dynamic 
exponent is z = 1.71(1), the S cLIHG cLS Z[ Qfi g to within uncertainty. This exponent 
is however considerably larger than the dynamic exponent z = 1.50(4) associated 
with finite-size scaling at the critical point. As discussed in 32J, the difference may 
reflect the existence of two time scales, one associated with the relaxation of the order 
parameter to its QS value, the other related to the finite-size lifetime. These two 
times scale in the same manner at absorbing phase transitions in models without a 
conserved density, such as the CP. These unexpected findings should motivate further 
study of the SRW and related models, with the goal of a complete and coherent picture 
of scaling in the CDP universality class. From the present vantage it appears that 
such a picture will be significantly more complex than for directed percolation. 
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